In [1]:
%pylab inline

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm

from scipy import stats
from itertools import product


Populating the interactive namespace from numpy and matplotlib
/Users/mbp_kirill/anaconda3/lib/python3.5/site-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.
  from pandas.core import datetools

In [2]:
wages = pd.read_csv('WAG_C_M.csv',';', index_col=['month'], parse_dates=['month'], dayfirst=True)

plt.figure(figsize(15,7))
wages.WAG_C_M.plot()
plt.ylabel('wages')
pylab.show()


свойства ряда:

  • повышающийся тренд
  • годовая сезонность
  • автокоррелированность
  • в начале ряда размах сезонных колебаний значительно меньше, чем в конце => дисперсия нестационарна

Можно с уверенностью утверждать, что ряд не обладает свойствами стационарности, но тем не менее воспользуемся критерием Дики-Фуллера, а также проведем STL-декомпозицию ряда:


In [3]:
plt.figure(figsize(15,10))
sm.tsa.seasonal_decompose(wages.WAG_C_M).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(wages.WAG_C_M)[1])


Критерий Дики-Фуллера: p=0.996132
<matplotlib.figure.Figure at 0x110682cc0>
  • тренд имеет простую, легко объяснимую инфляцией, структуру
  • сезонный профиль четко выражен

следующий шаг - стабилизация дисперсии:


In [4]:
wages['wages_box'], lmbda = stats.boxcox(wages.WAG_C_M)
plt.figure(figsize(15,7))
wages.wages_box.plot()
plt.ylabel('Transformed wages')
print("Оптимальный параметр преобразования Бокса-Кокса: %f" % lmbda)
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(wages.wages_box)[1])


Оптимальный параметр преобразования Бокса-Кокса: 0.270087
Критерий Дики-Фуллера: p=0.723160

И визуально, и воспользовашись критерием Дики-Фуллера, мы можем понять, что ряд по-прежнему не является стационарным. Но размах сезонных колебаний значительно уменьшился. Попробуем сезонное дифференцирование; сделаем на продифференцированном ряде STL-декомпозицию и проверим стационарность:


In [5]:
wages['wages_box_diff'] = wages.wages_box - wages.wages_box.shift(12)
plt.figure(figsize(15,10))
sm.tsa.seasonal_decompose(wages.wages_box_diff[12:]).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(wages.wages_box_diff[12:])[1])


Критерий Дики-Фуллера: p=0.008933
<matplotlib.figure.Figure at 0x1142fabe0>

Критерий Дики-Фуллера отвергает гипотезу нестационарности ряда на уровне значимости 0.05, но мы по-прежнему видим тренд. Попробуем добавить еще обычное дифференцирование:


In [6]:
wages['wages_box_diff2'] = wages.wages_box_diff - wages.wages_box_diff.shift(1)
plt.figure(figsize(15,10))
sm.tsa.seasonal_decompose(wages.wages_box_diff2[13:]).plot()   
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(wages.wages_box_diff2[13:])[1])


Критерий Дики-Фуллера: p=0.000000
<matplotlib.figure.Figure at 0x114230390>

Гипотеза нестационарности с уверенностью отвергается, и визуально ряд выглядит лучше. На трендовой компоненте не видно никакого систематического поведения - график колеблется в районе 0.

Перейдем к посмотрению модели. Для этого посмотрим на ACF и PACF полученного ряда:


In [7]:
plt.figure(figsize(15,8))
ax = plt.subplot(211)
sm.graphics.tsa.plot_acf(wages.wages_box_diff2[13:].values.squeeze(), lags=48, ax=ax)
pylab.show()
ax = plt.subplot(212)
sm.graphics.tsa.plot_pacf(wages.wages_box_diff2[13:].values.squeeze(), lags=48, ax=ax)
pylab.show()


  • максимальный сезонный лаг отличающийся от нуля отсутствует => Q = 0
  • последний несезонный лаг, при котором автокорреляция значима, 1 => q = 1
  • последний сезонный лаг, при котором частичная автокорреляция значима, 12 => P = 1
  • номер последнего несезонного лага, при котором частичная автокорреляция значима, 1 => q = 1 (значения остальных значимых лагов превышают сезонность, 12)

In [8]:
ps = range(0, 2)
d=1
qs = range(0, 2)
Ps = range(0, 2)
D=1
Qs = range(0, 1)

переберем все возможные значения параметров P, p, Q, q:


In [9]:
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)


Out[9]:
8

In [10]:
%%time
results = []
best_aic = float("inf")
warnings.filterwarnings('ignore')

for param in parameters_list:
    #try except нужен, потому что на некоторых наборах параметров модель не обучается
    try:
        model=sm.tsa.statespace.SARIMAX(wages.wages_box, order=(param[0], d, param[1]), 
                                        seasonal_order=(param[2], D, param[3], 12)).fit(disp=-1)
    #выводим параметры, на которых модель не обучается и переходим к следующему набору
    except ValueError:
        print('wrong parameters:', param)
        continue
    aic = model.aic
    #сохраняем лучшую модель, aic, параметры
    if aic < best_aic:
        best_model = model
        best_aic = aic
        best_param = param
    results.append([param, model.aic])
    
warnings.filterwarnings('default')


wrong parameters: (0, 0, 0, 0)
CPU times: user 2.5 s, sys: 72 ms, total: 2.58 s
Wall time: 1.43 s

In [11]:
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
print(result_table.sort_values(by = 'aic', ascending=True).head())


     parameters        aic
4  (1, 0, 1, 0)  38.903410
2  (0, 1, 1, 0)  39.603519
6  (1, 1, 1, 0)  40.900114
3  (1, 0, 0, 0)  41.921386
1  (0, 1, 0, 0)  42.608284

Лучшая модель (минимальный aic = 38.9, parameters = (1, 0, 1, 0)):


In [12]:
print(best_model.summary())


                                 Statespace Model Results                                 
==========================================================================================
Dep. Variable:                          wages_box   No. Observations:                  294
Model:             SARIMAX(1, 1, 0)x(1, 1, 0, 12)   Log Likelihood                 -16.452
Date:                            Thu, 17 Aug 2017   AIC                             38.903
Time:                                    21:04:36   BIC                             49.954
Sample:                                01-01-1993   HQIC                            43.329
                                     - 06-01-2017                                         
Covariance Type:                              opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.1745      0.045     -3.867      0.000      -0.263      -0.086
ar.S.L12      -0.1405      0.046     -3.033      0.002      -0.231      -0.050
sigma2         0.0658      0.004     15.140      0.000       0.057       0.074
===================================================================================
Ljung-Box (Q):                       54.52   Jarque-Bera (JB):                46.95
Prob(Q):                              0.06   Prob(JB):                         0.00
Heteroskedasticity (H):               1.56   Skew:                             0.16
Prob(H) (two-sided):                  0.03   Kurtosis:                         4.98
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

Остатки:


In [13]:
plt.figure(figsize(15,8))
plt.subplot(211)
best_model.resid[13:].plot()
plt.ylabel(u'Residuals')

ax = plt.subplot(212)
sm.graphics.tsa.plot_acf(best_model.resid[13:].values.squeeze(), lags=48, ax=ax)

print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(best_model.resid[13:], 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(best_model.resid[13:])[1])


Критерий Стьюдента: p=0.137126
Критерий Дики-Фуллера: p=0.000080

Остатки несмещены (подтверждается критерием Стьюдента) стационарны (подтверждается критерием Дики-Фуллера и визуально), неавтокоррелированы (подтверждается критерием Льюнга-Бокса и коррелограммой).

Посмотрим, насколько хорошо модель описывает данные:


In [14]:
def invboxcox(y,lmbda):
   if lmbda == 0:
      return(np.exp(y))
   else:
      return(np.exp(np.log(lmbda*y+1)/lmbda))

wages['model'] = invboxcox(best_model.fittedvalues, lmbda)
plt.figure(figsize(15,7))
wages.WAG_C_M.plot()
wages.model[13:].plot(color='r')
plt.ylabel('wages')
pylab.show()


Прогноз:


In [15]:
wages2 = wages[['WAG_C_M']]
date_list = [datetime.datetime.strptime("2017-07-01", "%Y-%m-%d") + relativedelta(months=x) for x in range(0,24)]
future = pd.DataFrame(index=date_list, columns= wages2.columns)
wages2 = pd.concat([wages2, future])
wages2['forecast'] = invboxcox(best_model.predict(start=294, end=317), lmbda)

plt.figure(figsize(15,7))
wages2.WAG_C_M.plot()
wages2.forecast.plot(color='r')
plt.ylabel('wages')
pylab.show()


Прогноз выглядит адекватным, соответствует тренду и передает информацию о сезонности


In [ ]: